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A low frequency stochastic background of gravitational waves may be detected by pulsar timing 
experiments in the next five to ten years. Using methods developed to analyze interferometric 
gravitational wave data, in this paper we lay out the optimal techniques to detect a background of 
gravitational waves using a pulsar timing array. We show that for pulsar distances and gravitational 
wave frequencies typical of pulsar timing experiments, neglecting the effect of the metric perturbation 
at the pulsar does not result in a significant deviation from optimality. We discuss methods for 
setting upper limits using the optimal statistic, show how to construct skymaps using the pulsar 
timing array, and consider several issues associated with realistic analysis of pulsar timing data. 

PACS numbers: 



I. INTRODUCTION 

The search for gravitational waves is at the forefront 
of current fundamental physics research. The direct de- 
tection of gravitational waves will usher in a new era in 
astronomy and astrophysics. Gravitational waves will re- 
veal information about black holes, supernovae, and neu- 
tron stars that cannot be gleaned from electromagnetic 
observations. Furthermore, the detection of a gravita- 
tional wave background will open an observational win- 
dow onto a time in the early universe before recombina- 
tion, prior to which the universe is opaque to electromag- 
netic waves. The scientific rewards for such a detection 
would be truly exceptional. Several international efforts 
are underway to detect gravitational waves and two of 
these efforts are expected to result in the detection of 
gravitational waves in the next 5 to 10 years: Interfer- 
ometric ground-based gravitational wave detectors and 
pulsar timing observations. 

Neutron stars radiate powerful beams of radio waves 
from their magnetic poles. If a neutron star's magnetic 
poles are not aligned with its rotational axis, the beams 
sweep through space like the beacon on a lighthouse and 
the neutron star is said to be a pulsar. If the Earth lies 
within the sweep of a pulsar's beams, the star is observed 
as a point source in space emitting short, rapid bursts of 
radio waves [Ij. Due to their enormous mass, neutron 
stars have a very large moment of inertia and the radio 
pulses we observe arrive at a very constant rate. Pul- 
sar timing experiments exploit this regularity [31 HI E] . 
Fluctuations in the time of arrival of radio pulses, af- 
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ter all known effects have been subtracted out, could be 
due to the presence of gravitational waves. Since the 
1970s, when these ideas were first conceived, pulsar tim- 
ing precision has improved dramatically. Several known 
pulsars can now be timed with a precision of about 1 
micro-second, and a handful can be timed with a pre- 
cision around several hundred nanoseconds Recent 
work [7 has shown that the presence of nanohertz gravi- 
tational waves could be detected by observing 20 pulsars 
with timing precisions of 100 nanoseconds over a period 
of 5 to 10 years. Non-detection would still improve cur- 
rent bounds on the low frequency stochastic gravitational 
wave background [8]. 

Gravitational waves from supermassive black hole bi- 
nary systems could be detected via pulsar timing obser- 
vations [HI Uni im eh CSl- in addition, pulsar timing 
has the potential to measure the polarization proper- 
ties of gravitational waves which could confirm (or even 
change) the current theory of gravity \^ [T5]. Gravita- 
tional wave observations in the nanohertz band could also 
yield information about the early universe |16j . Cosmic 
strings, line-like topological defects, could produce grav- 
itational waves in the nanohertz band. Cosmic strings 
can form during phase transitions in the early universe 
due to the rapid cooling that takes place after the Big 
Bang [T71 [T5J [12]. Cosmic string production is generic 
in supersymmetric grand unified theories [21)]. Addi- 
tionally, in string theory motivated cosmological mod- 
els cosmic strings may also form (dubbed cosmic super- 
strings to differentiate them from field theoretic cosmic 
strings) US HU US [Ml I25J US]- Cosmic strings and 
superstrings are expected to produce a background of 
stochastic gravitational waves and bursts of gravitational 
waves [27, 28 , 29, 30, 31, 32 that could be detected using 
pulsar timing observations. Pulsar timing observations 
are already producing some of the most interesting con- 
straints on cosmic string models and a detection would 
have profound implications [5111^ . 

Recently Lommen |33j produced an upper limit on the 
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stochastic gravitational wave background using observa- 
tions of three miUisecond pulsars spanning 17 years. The 
methods used by Lommen were based on those devel- 
oped by Kaspi and collaborators "31] and have been the 
subject of some criticism in the literature [8j i30j. More 
recently Jenet and collaborators [H |8] developed a new 
technique for gravitational wave stochastic background 
searches in pulsar timing data and applied it to Parkes 
Pulsar Timing Array data [351 ES] • 

In this paper we consider optimal strategies for extrac- 
tion of a gravitational wave stochastic background signal 
using data from a pulsar timing array. Our methods are 
based on those developed for and used in ground based in- 
terferometric gravitational-wave detectors such as LIGO 
and Virgo ^7] [SSJ [Ml Sni IHl iU il Hi Hi], and im- 
prove on existing methods in several ways. In Section [ll] 
we write the redshift of pulsar signals induced by pass- 
ing gravitational waves, first derived by Detweiler [3], 
in a coordinate-independent way more suitable for our 
analysis, and discuss its form in the frequency domain 
including the long- wavelength limit. In Section |III| we 
construct the optimal cross-correlation filter for a pulsar 
pair by maximizing the signal to noise. We find that 
the overlap reduction function is well approximated by 
a constant, or equivalently, that the metric perturbation 
at the pulsar can be neglected for values of pulsar dis- 
tances and gravitational wave frequencies typical of pul- 
sar timing experiments, without significant losses in sen- 
sitivity. In Section jlll] we also show how to construct the 
optimal combination of cross-correlations of pulsar pairs 
in a pulsar timing array and include a more sophisti- 
cated derivation of the optimal detection statistic based 
on the likelihood ratio. In Section |IV| we discuss up- 
per limit and detection methods. In Section |V] we show 
how to construct skymaps using pulsar timing data — 
a pulsar timing radiometer. In Section |VI| we discuss 
several important issues relating to the realistic analysis 
of pulsar timing data, including the Lomb-Scargle peri- 
odogram for power spectrum estimation of unevenly sam- 
pled data, and optimal procedures for computing Fourier 
transforms. We conclude in Section VIII. Lommen, Ro- 
mano and Woan [IS] will extend our work using a likeli- 
hood based approach developed in [ITj, and consider the 
case of stochastic backgrounds that are loud compared 
to the noise, closely examine time-domain implementa- 
tions of the optimal statistic, and provide a detailed com- 
parison of the optimal statistic described here with the 
methods of Jenet and collaborators [ZIIH]- 



II. THE SIGNAL 

Gravitational waves affect pulsar timing measurements 
by creating perturbations in the null geodesies that the 
radio signals emitted from the pulsar travel on j3]. In 
this section we will describe the relationship between the 
metric perturbation and the signal measured in pulsar 
timing experiments. 



A metric perturbation in a spatial, transverse, traceless 
gauge has a plane wave expansion given by |40j 

where / is the frequency of the gravitational waves, k — 
27r/f2 is the wave vector, f2 is a unit vector that points 
along the direction of travel of the waves, i,j — x,y,z 
are spatial indices, and the index A = +,x labels polar- 
izations. The polarization tensors efj{fl) are 

efJfl) = rhiirij — •hi'hj, 



(2a) 
(2b) 



where 



Cl = (sin cos 0, sin 6 sin cos 9) , (3a) 
m (sin 0, — cos 0, 0) , (3b) 
h — (cos 9 cos cos 9 sin </i, — sin 9) . (3c) 

Now consider the metric perturbation from a single 
gravitational wave traveling along the z-axis so that f2 — 
z. The metric perturbation is given explicitly by 

/oo 

= h,j{t-z). (4) 
The physical metric due to the perturbation is given by 



gab = r/ab + habit-z) 



-10 

1 + /1+ /ix 

/ix 1-/1+0 

1 



, (5) 



where rjab = diag{— 1, 1, 1, 1} is the Minkowski metric, 
a, b are spacetime indices, and 



i+^x{t - z) 



d/e^2^/(*-^)/i+,x(/,5). (6) 



In this background, a pulsar emitting pulses at frequency 
vq and direction cosines a, [3, and 7, with respect to 
the X-, y-, and z-axes, respectively, will be observed to 
change its frequency in the solar system reference frame 
according to [3] 



z{t,z) 



o?- - 3"^ ad 



where h'^ x ; ^+ x the gravitational wave strains at 
the solar-system barycenter and the pulsar, respectively. 
This central result was obtained by Detweiler |3], who 
generalized a result of Estabrook and Wahlquist [5] to 
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include both gravitational wave polarizations and for pul- 
sars at arbitrary locations. They in turn based their cal- 
culation on an earlier one by Kaufmann [15] . A detailed 
derivation of this result is provided in Appendix \K\ 

Looking at Eq. ^ (and as shown in Appendix |Bp we 
can write the redshift z(t, Cl) of signals from a pulsar in 
the direction of the unit vector p produced by a gravita- 
tional wave coming from the direction as, 



2i + n-p 



-Ah. 
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where 



Ah, 



(8) 



(9) 



is the difference in the metric perturbation traveling 
along the direction Cl at the pulsar and at the center 
of the solar system. The vectors {te,Xe) and {tp,Xp) give 
the spacetime coordinates of the solar-system barycenter 
and the pulsar, respectively. The metric perturbation at 
each location takes the form, 



hij{t, fl) = 



A 



dfe''^f^'-''-'°'>hA{f,Cl)ef^{Cl), 



(10) 
for a fixed Xq. 

We choose a particular coordinate system by placing 
the solar-system barycenter at the origin and the pulsar 
some distance L away. With these conventions 



tc- 

0, 

Lp. 



L = t- 



(11) 
(12) 
(13) 



If assume that the amplitude of the metric perturbation 
is the same at the solar-system barycenter and the pul- 
sar then we can use Eq. (10) to write out Ahij in our 
coordinate system as 



Ah, 



-i27r/L(l+0-j5) 



1 



xyhAimefm) 



A 



Ah,j{t,n). 



(14) 



Ultimately, we will be interested in the Fourier transform 
of this quantity which is simply 



Ah„{f,n) = e 



-j:27r/L(l+f2-p) _ 



A 



hAif,n)efjn) 



(15) 



We can then write the Fourier transform of Eq. n8n as 



~i2TTfL{l+n-p) 



A 



hAif,n)F^{n), 



(16) 



where we have defined 



1 



p p' 



2l + fl-p 



(17) 



As shown in Appendix |B] the total redshift is given by 
summing over the contributions coming from gravita- 
tional waves in every direction: 



dni{f,n), 



(18) 



and similarly for z(t). 

In fact, it is not the redshift, but a related quantity 
called the residual that gets reported in pulsar timing 
measurements. The residual, R{t), is defined as the inte- 
gral of the redshift: 



m 



dt' z{t'). 



(19) 



This simple relationship gives us the freedom to develop 
the data analysis for either variable and we henceforth 
limit our attention to the redshift, but the results here 
can be phrased in terms of the residual with minimal 
effort. 

In the literature, searches for gravitational waves using 
pulsar timing data are typically performed in the time do- 
main. The (unknown) metric perturbation at the pulsar 
in, say, Eqs. (|7| or ^ is neglected because one can treat 
it as another noise term which averages to zero when 
performing correlations between measurements of differ- 
ent pulsars. In the frequency domain this is unnecessary. 



Eq. (16) does not depend explicitly on the metric per- 



turbation at the pulsar, rather the dependence is all in 
a distance and frequency dependent phase factor. It is 
then conceivable that if we could determine the distance 
to a pulsar L with sufficient accuracy we could use the 
metric perturbation at the pulsar to improve the sensi- 
tivity of our searches. Unfortunately, such measurements 
of pulsar distances are unavailable. We will show in Sec- 



tion III A however, that for the case of a stochastic back- 



ground search in pulsar timing data the phase factor can 
be neglected without any significant loss in sensitivity. It 
is unclear whether this is true for other types of gravita- 
tional wave searches. 

From the ground-based interferometer perspective 
Eqs. or (jsj) are somewhat counter-intuitive. This diffi- 
culty arises from the factor oi\ + U,-pm the denominator; 
in Appendices]^ and [B] we show explicitly how this factor 
enters the expression. When Cl-p = ±1, i.e. the grav- 
itational wave and the pulsar directions are parallel or 



anti-parallel, Eqs. (16) and (17) lead to no redshifting of 



the pulsar signal for completely different reasons. When 
they are parallel the reason is the transverse nature of 
gravitational waves, and when they are anti-parallel it is 
because the pulsar signals "surf" the gravitational waves. 
Our surprise is a result of our long-wavelength limit in- 
tuition. 



4 



Eq. ( 16 ) has an obvious long- wavelength limit. We can 



use this limit to compare the form of our results with 
those of ground-based interferometers such as LIGO. 
When 27r/i ^ 1 we can Taylor expand the exponential 
and to first order Eq. becomes 

(20) 



~z{f,n) « -z^/WX^-4(/, j^)e^J(f^). 



Typical values of / are in the range 1/10 yr^^ to 10 yr^^. 
Typical values of the Earth pulsar distance L are in the 
range 100 ly to 10"*^ ly. This means fL is in the range 
10 to 10^ and pulsar timing experiments are never in the 
long wavelength limit. However, the Taylor expansion 
can also be done for large fL when the angle between Ct 
and p is sufficiently close to tt. In this case the pulsar sig- 
nals can "surf" the gravitational waves and not undergo 
redshifting. Writing that angle as tt — e with e ^ 1 , then 
the Taylor expansion is also valid when e ^ (7r/i)~^/^. 
Taking the inverse Fourier transform of Eq. ( pO| yields 

L 



m 



-fp'hij{t,Xc), 



(21) 



which is the projection of the time derivative of the met- 
ric perturbation at the solar-system barycenter onto the 
unit vector that points to the pulsar. Note that unlike 
Eq. ([s]) , this equation no longer depends on the direction 
of the gravitational wave and can be expressed in terms 
of the full metric perturbation (derivative). For the case 
of ground based interferometers the signal, the so-called 
strain, is proportional to the difference in length of the 
two arms because the signal at the dark port of the inter- 
ferometer depends on that difference. If the arms point 
in the directions of the unit vectors X and Y the strain 
is given by 



h{t) 



(22) 



which is the metric perturbation hij{t, x) projected onto 
the difference of the arms. 



III. DETECTION STATISTIC 

With an understanding of the signal in hand we now 
turn our attention to developing an optimal detection 
strategy. In this section we will first derive the optimal 
cross-correlation statistic for a single pulsar pair using ar- 
guments based on maximizing signal to noise ratio. We 
will then determine the best way to combine measure- 
ments from multiple pulsar pairs to obtain the most con- 
straining upper limit. This section will conclude with 
a more sophisticated derivation of the optimal detection 
statistic based on the likelihood ratio. 



A. The optimal filter 

In this section we will derive the optimal filter for de- 
tecting a stochastic background of gravitational waves 



from the cross-correlation of redshift measurements of 
two different pulsars. This problem was addressed in 
detail by Allen and Romano [ID], for the case of inter- 
ferometric gravitational wave detectors and our analysis 
follows theirs closely. 

Consider the signals from two pulsars 



si(i) = zi(i)+ni(t), 

32{t) = Z2{t)+n2{t), 



(23) 
(24) 



where Zi(t) is the redshift and ni{t) is the noise intrinsic 
in the measurement. Throughout this work we will as- 
sume that each ni(t) is stationary and Gaussian, and is 
greater in magnitude than the redshift. Additionally we 
assume that 



{z^{t)) 

(ni(i)n2(t)> 
{n,{t)z,{t)) 



0, 
0, 
0, 
0, 



(25) 



for all i and j, where the angle brackets denote an expec- 
tation value. 

A stochastic background will show up in the data as 
correlated noise between measurements with different de- 
tectors. Our goal is to find a filter, Q{t — t'), that opti- 
mizes the cross-correlation statistic 



S = 



T/2 



dt 



T/2 



T/2 



T/2 



dt' si{t)s2{t')Q{t-t'), (26) 



where T is the observation time. We will define the op- 
timal filter to be the Q{t — t') that maximizes the signal 
to noise ratio (SNR) 



SNR = 



(27) 



where /i and a are the mean and square root of the vari- 
ance, respectively, associated with the cross-correlation 



signal defined in Eq. ( 26 1 . 

We start by assuming that the observation time is 
much greater than the separation of the two detectors 
and extend the limits of the integral over dt' to ±00. 
Technically our assumption is not correct because pul- 
sars are typically separated by distances far greater than 
the observation time. Later we will see that neglecting 
the phase terms that correspond to the metric perturba- 
tion at the pulsar location is an excellent approximation. 
In effect this makes our detectors co-located though not 
co-aligned, and our assumption about the observation 
time is appropriate. We work in the frequency domain 
so that Eq. ( 26 1 becomes 



drsT{f-fr4{f)h{f)Q{f), (28) 



where (5t(/ ~ /') is the finite-time approximation to the 
delta function 



Srif) 



sin(7r/r) 



(29) 
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The mean of the cross-correlation is 



M ^ (5) = 



df 



df 5T{f ^ f){ziu)uf))Q{n- 

_ J 

(30) 

Because of Eqs. p^, (jlSl), and (pSl), taking the 



expectation value above requires us to evaluate 
{h\{f, h)hA'{f' , ^'))- The assumptions that the stochas- 
tic background is stationary, unpolarized and isotropic 
lead us to take 

{h\{f, n)hA'{f, n')) = 6\n, Ci')5aa'5U ~ DhU), 

(31) 

where H{f) = H{—f) is the gravitational wave spec- 
trum. H{f) is related to rigw(/) through 



^gw(/) — 

where pcrit — Stt/SHq and 



1 dp. 



gw 



Pcrit rfln/' 



1 

32^ 



(habit, X)h-''{t,x)), 



(32) 



(33) 



is the energy density in gravitational waves. It fol- 
lows from the plane wave expansion Eq. ([T]) along with 
Eqs. pT) and ([32| in Eq. (|33l that 



^(/) = ^i/r'f^gw(i/i), 



(34) 



and therefore 



3F, 



327r3 

xi/r'f^gw(i/i), 



(35) 




fL > 10 



= 0' ^ - 0.86 Pulsar timing experiments 




r = -0.29, ^ = It/2 



FIG. 1; Plot of the full overlap reduction, Eq. (381, along 



with the approximation Eq. ( 39 1 for two pulsars a distance 
L from the solar-system barycenter. The overlap reduction 
function is a function of fL. The top two (blue) curves show 
Eq. fMl with /? = 3/47r (solid line) and Eq. fMl (dashed 



line) for two pulsars at an angle ^ — tv/8 as a function of fL. 
The middle (red) and bottom (green) curves show the same 
quantities for two pulsars at ^ « 0.86 and ^ = 7r/2 respec- 
tively. The smallest value of the frequency /min ~ 0.1 yr"'^ 
and the closest pulsars used in timing experiments are at a 
distance of Lmin ~ 100 ly so that fL > 10. This range of fL 
puts pulsar timing experiments in the regime where Eq. 1 39 1 
is an excellent approximation to Eq. (381 and we are justi- 



fied in throwing out the pulsar term while remaining close to 
optimal. 



which is sometimes written in terms of the characteristic 
strain 



27r2 /2"s* 

The expectation value we set out to evaluate is then 



(36) 



3gg 1 
327r3 f3 



s{f-f)\f\-'n,u\f\n\f\), 



where we defined 

r(l/l)= /? E 



(37) 



(38) 

the pulsar timing analogue of the overlap reduction func- 
tion [OT, which has a normalization factor /3. The nor- 
malization is chosen so that r(|/|) = 1 for coincident, 
co-aligned detectors. As we show below, pulsar timing 
experiments are in a regime where the exponential factors 



in Eq. (38 1 can be neglected. In this situation, which we 



will assume henceforth, the normalization factor is easy 
to determine and we have that 



To 



= 3 



-T 

A 

1 

3 + 



dVLF('{p)F^{pL) 
1-cose 



1 — COS ^ 



(39) 



where ^ = cos~-^(pi ■ ■p-i) is the angle between the two 
pulsars. This quantity is proportional to the Hellings and 
Downs curve [49]. A detailed derivation of this result is 
provided for completeness in Appendix |C 1[ 

The rationale given in the literature for throwing out 
the pulsar term in Eq. ([7|, or equivalently Eqs. ([8| and 
([9]), is that the unknown metric perturbation at the pul- 
sars can be thought of as a kind of noise term which aver- 
ages to zero when performing a correlation between dif- 
ferent pulsars. The equivalent procedure in the frequency 
domain is to neglect the phase factors in Eq. (16), or in 
terms of our optimal filter, approximating Eq. (38 1 with 
Eq. ( 39 ) . The regime where the approximate Eq. ( 39 ) is 
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valid, is helpful in quantifying the accuracy of the ratio- 
nale. Figure [T| shows the overlap reduction function for 
two pulsars a distance L from the solar-system barycen- 
ter. Since the distance to both pulsars is the same, the 
overlap reduction function can be written as just a func- 
tion of JL. The top two (blue) curves show Eq. ( 38 ) with 
(3 = 3/47r (solid line) and Eq. (39 1 (dashed line) for two 
pulsars at an angle ^ = 7r/8 as a function of fL. The 
middle (red) and bottom (green) curves show the same 
quantities for two pulsars at ^ ~ 0.86 and ^ = 7r/2 re- 
spectively. As discussed in the last section the smallest 
value of the frequency /min ~ 0.1 yr~^ and the closest 
pulsars used in timing experiments are at a distance of 
imin ~ 100 ly so that fL > 10. As shown in Fig. [T|this 
range of fL puts pulsar timing experiments in the regime 
where Eq. ( 39 1 is an excellent approximation to Eq. (|38| , 



and we can neglect the pulsar term while remaining close 
to optimal. 

Returning to Eq. ( 30 ) , we now have 



3gg 
327r3 (3 



dfifr^n^uimm 

-\ttJ df\f\-^n^^m- (40) 



With the assumption that the noise is much greater 
than the signal, the variance, cr^, depends only on the 
statistical properties of the noise in each detector. We 
have 



T 
4 



dfP,{\f\)P2{\f\) Q{f) 



where 



(nUfW)) = ^S{f - f')mf\), 



(41) 



(42) 



is the (one-sided) noise power spectrum. 

With ^ and in hand, we next define a positive- 
definite inner product using the noise power spectra of 
the two detectors 

/oo 
dfA*{f)Bif)P,i\f\)P2{\f\)- (43) 
-OO 



With this definition it is easy to see that 



T 
4 



»gw(|/|)ro 

|/PPl(|/|)P2(|/|) 



(q,q) , 



(44) 
(45) 



from which it follows from the definition of SNR in 



Eq. (|27| and Schwartz's inequality that the optimal filter 

f^gw(|/|)ro 



is given by 



Qif) = X 



\f\'Pi{\f\)P2m" 



(46) 



for some normalization constant, x- Our primary interest 
will be in stochastic backgrounds with power law spec- 
tra, rigw(/) = i^a/" (for constant fia). In that case the 
normalization constant for the optimal filter, Qa{f), is 
chosen so that 



(47) 



where Tq is some arbitrary constant with dimensions of 
time. From Eq. ( 44 1 it follows that 



Y — \ ln ?r 



df 



W(|/|)P2(|/|) 



Finally, we can compute 



SNR « Mri/2 

47r2 



df 



»^w(l/l)rg 



1/2 



(48) 



(49) 



The differences between these results and those for in- 
terferometers can all be traced to the differing overlap 
reduction function r(/) « Fq. The normalization of Fq 
means that the maximal SNR (for co-incident, co- aligned 
detectors) is only 5 /6 of that obtainable from interferom- 
eters, assuming the noise power spectra are the same in 
each case. 



To construct the optimal filter, Eq. (46 1, the noise 
power spectra for the two pulsars Pi (|/]) and -P2(|/|) 
must be determined. These can either be modeled, or 
measured with the methods described in Section IVII 
Once constructed the optimal filter can be applied in 
the frequency domain. Section VI gives a prescription 



for taking Fourier transforms of unevenly sampled data. 
The optimal filter can also be inverse Fourier transformed 
and the correlation performed in the time domain. It is 
unclear which of these two methods is more robust and 
the authors of j46l will explore the time-domain approach 
in detail. 



B. The pulsar timing array 

The question we would like to address in this section 
is: Given redshift measurements from N different pulsars 
(which each have a different noise profile), what is the 
best way to combine those measurements to produce the 
most constraining upper limit? One can consider the 
cross-correlations between any even number of detectors, 
but it has been shown [40l [50] that the optimal choice 
is the combination of pairwise cross-correlations. As it 
turns out, the solution to this problem also solves the 
problem of non-stationarity in the noise power spectra 
over periods longer than the typical observation time, T. 

First let 



(50) 



be riij measurements of the cross-correlation between 
the i-th and j-th pulsar. We will assume that each 
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measurement is taken with an optimal filter normalized 
so that while searching for a background of the form 



for some constants, which has mean 



(51) 



and variance 



where Tq is an arbitrary constant introduced for dimen- 
sional reasons. Each measurement therefore has the form 



Mm = (A) = A* 



(60) 



I I 



2 2 

ij ij 



df 



(52) 




(61) 



with 



(ij) 



Qkif) 



(ij) 



Xk 



»gw(|/|)fa"'ro 

\f\'PrA\f\)P,A\f\)' 



(53) 



where ^^^^Fq is the overlap reduction function of the (ij) 
pulsar pair, Si^k is the k-th measurement of the signal 
from the i-th pulsar, and Pi^k is the associated noise 
power spectrum. Additionally, 



df 



pp.AifmAifi) 



(54) 

where ^"^-"Tk is the observation time of the fc-th measure- 
ment of the (ij) pulsar pair. Our task is to combine the 
^*^-'S'fe in a way that optimizes SNR. The first step is to 
form the sample mean for each set of measurements 



n - ■ ^ — ^ 



(55) 



k=l 



which is both an unbiased estimator and random vari- 
able. It therefore has a mean 



and a variance 



t2. = /(■'J)r,2 



(56) 



(57) 



where 



"21 (ij)^! /-oo 

fe=i 



df^^^hi 



pp^A\f\)PoA\f\y 



(58) 

The next step is combine the sample mean for each set of 
measurements into a single estimator we can use to de- 
termine an upper bound on Qq, and hence figw = ^af"'- 
We do so by introducing an unbiased estimator consisting 
of a weighted average of the sample means 



1=1 j<i 



I I 



(59) 



EE A. 

i—1 j<i 



The object is to now determine the A.^- that maximizes 
the SNR of /t. The (squared) SNR of jl is 



SNR2 = ^2 




EE 4- 

1=1 j<i 



(62) 



To find the A^j that maximize the SNR, we exploit the 
same trick that led us to the optimal filter. Namely, we 
introduce an inner product 



( / 

i—1 j<i 

which allows us to write 

from which it follows that choosing 



(63) 



(64) 



(65) 



maximizes the SNR. The optimal statistic, choosing 



A, 



fjjj , is then given by 



; I 



2—1 j<i 



-'opt 



EE'^.; 

i—1 j<i 



i=l j<i 



k=l 



I I 

EE 

i—1 j<i 



(66) 



Because the estimator defined in Eq. ( 59 ) is unbiased and 
defined so that /i = (Sopt) — ^aTo, the estimate of Cta is 
found using 



So 



pt 



(67) 
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where Sopt is the measured value of the optimal statistic. 
The expected variance of S'opt follows from Eq. (61 ), 



Qk but cancels in Eq. ( 67 1 so that the point estimate 



of JIq is independent of Tq. 



/ I 

EE' 

i—l j<i 



(68) 



Aside from maximizing the SNR, the linear combina- 
tion of sample means that forms the optimal statistic 
in Eq. (66 1 serves two important and related purposes. 



First of all, as mentioned at the beginning of this sec- 
tion, weighing each (*^)/t by the inverse of the squared 
variance means that less noisy measurements (those with 
smaller variances) contribute more to the sum, which 
helps minimize the effect of long term non-stationarity. 
This is augmented by the normalization convention we 
chose in Eq. (511 for the mean of each measurement. Us- 
ing Eg. ( 45 1 and Eq. ( 46 ) with the A that follows from 



Eq. (51 1, we see that ^ cx (^-^^T, and so longer ob- 



servation times are also favored in the sum. 



C. Computational procedure 

In this subsection we describe how the quantities nec- 
essary for a stochastic background search are computed. 
The goal is to produce ameasurement of the optimal 
statistic, S'opt, using Eq. (66). The optimal statistic can 



then be used to make detection or upper limit statements 



(see Section IV I 



First the power spectra spectra for each pulsar (and 
each stretch), Pi^k{\f\), must be determined. The spec- 
tra can either be modeled or measured with the methods 



described in Section VI Then the overlap reduction func- 
tions, '•^-'^Fo, need to be computed for each pulsar pair. 
To optimize the statistic for particular spectra the value 
of a (in rigw(/) = f^a/") needs to be chosen. The nor- 
malizations, ^*-'^Xfc, can then be computed using Eq. (54 1. 



The normalizations a llow us to compute the variances. 



a , given by Eq. (58 1, in the numerator and dcnomi- 



nator of Eq. (66 1, as well as the filters, ^*-''Qfc(/), through 
Eq. (53 1. Note that the unknown factors of ^la cancel 



everywhere: From Eq. (54 1 it is easy to see the normal- 
ization '■'-'-'xfe oc r^Q^, so there is a cancellation a factor 
of in Eq. (Issf, and a factor of fll in Eq. With 



these quantities in hand the cross-correlations, '^'^^'Sk, in 
Eq. ( 52 1 can be computed by taking Fourier transforms 



of the data (see Section VI I . Alternatively, a set of time- 
domain filters, ^^^^Qk{t), can be created by taking inverse 



Fourier transforms of Eq. ( 53 1 and applied to the data in 



the time domain using Eq. ( 26 1 . 

Note that there is no dependence on the arbitrary 



constant Tq introduced in Eq. (51 1 for dimensional rea- 
sons. The ^'-'-'xfc Eire linear in Tq and enter the variances 
quadratically (see Eq. (58 1). The dependence cancels in 
Eq. ( 66 1 because it is present in both numerator and de- 



D. Likelihood approach 

The detection statistic that has been derived is also 
an optimal statistic in the sense that it is the logarithm 
of the likelihood ratio, at least in the limit where the 
expected signal is smaller than the noise, and therefore 
it is the optimal statistic in both the Bayesian sense and 
by the Neyman-Pearson criterion. This section is based 
on the likelihood analysis of |17], generalized to consider 
multiple detector pairs. 

As we did previously, we assume that the noise is sta- 
tionary and Gaussian, as is the stochastic background. 
For any given pulsar i we assume that there are discrete 
samples of data which forms a vector s^. Although the 
discussion below does not place requirements on the data 
sampling, we will assume that the observations of the pul- 
sars all involve the same number of points N at the same 
evenly spaced sampling interval so that sample j of pul- 
sar i is Si[j] = Si(jAt) where At is the sampling interval. 
This signal vector is the sum of a noise vector and the 
redshift vector z^, Sj ~ Zi+Ui. The data is a combination 
of two random processes: the instrumental noise and the 
contribution from the stochastic background. The auto- 
correlation matrix = (s| Sj) is an x matrix 
which contains both of these contributions and, since we 
assume Gaussian noise and stochastic background, this 
matrix completely characterizes the distribution of the 
data. As we did previously, we assume that the measure- 
ment noise in a pulsar observation is independent of the 
noise in the observations of other pulsars; the stochastic 
background, however, is correlated amongst the pulsar 
signals. This correlation is characterized by the stochas- 
tic background correlation matrix e^Sy — {zl^Zj). Here 
e is an order parameter which we will use to expand the 
probability distribution in powers of the small stochastic 
background signal. It can also be interpreted as an over- 
all amplitude parameter of the stochastic background. 
The probability distribution for the collection of all pul- 
sar observations is given by a multidimensional Gaussian 
distribution 



p(x|e) 



i/det(27rS) 



where 



exp I 



Sl 
S2 



2^ 



(69) 



(70) 



nominator. Tq also enters S'opt linearly through '•^-'^Xfe in is a column vector formed from all of the data vectors 
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and 





Ri e^Si2 • 








£^821 R2 






s = 




ni 


(71) 



is the correlation matrix for the collective observation 
vector X. In this weak signal limit we find 



R]"^ 
R:7^ 



Rr^ 



R]^ ^"812" R2 ^ 

R2 ^ • S2I • R;^ ^ 

. R; ^ ■ ■ Ri ^ R; ^ ■ ^'2 ■ R2 ^ 



Rl ^ • Si; • R; ^ 

Ro • S2; • R; 



Rl ^ • Sito • Rm^ • Smi • Rj ^ Rj ^ • Sim • R^^ ' Sm2 ' R2 ^ 

m — 1 m=l 

mT^l TO7tl,2 

; i 

R2 ^ ■ S2m ■ Rm^ ■ Sml " Rl ^ R2 ^ " " Rm^ ' Sm2 ' R2 



Rl ^ • Sim • Rto^ ■ Sto; • R; ^ 

I 

R2 ^ ■ S2m • Rto^ • Sto; • R; ^ 



m=l 



m=l 
m^2 



m=l 
m#2,; 



^ ^ R; ^ ' Sim • Rm^ ■ Sml ' Ri ^ ^ ^ R; ^ ' • R^^ • 8^2 ' R2 ^ ' ' ' ^ ^ R; ^ ' 8;^ • R^"'^ • 8^; ' R; 



m—l 
m^l,l 



m—l 



m—l 
m^l 



and 



(72) 



In det S = ^ In det Rj + ^ ^ tr(Rri . g . j^-i . g^. _^ ^^^e^ 

i=l i=l j<i 



(73) 



The logarithm of the likelihood ratio is 

InA = lnp(x|e) - lnp(x|0) 

I I 



i=l j<i 



1 ' 



i=l 



; I 



J2 tr(Rr' • Si,- • R7^ • 8,,) - 2 ^ ^ JR (si • R71 . 8iTO • R-^ • S^,- • Kj' ■ s,) 

j<i i<i m=l 



+0(e«) 



= e25-ie4AA2 + 0(e6). 



(74) 



This is the optimal detection statistic for a weak stochas- 
tic background. We have identified <S as the O(e^) term 
and — 2A/'^ as the 0{e*) term of the log-likelihood ratio. 

The locally optimal detection statistic is obtained in 
the e limit; it is the leading 0{e'^) term: 



, InA 

hm ^5- 



= '5 = EE^(«I-R-»"'-s^^-i^7'-«^) 



i=l j<i 



(75) 



Although this presentation has been described in terms 
of observational vectors in the time-domain, the deriva- 
tion of the likelihood ratio has not explicitly required 
this choice of basis. It is convenient to perform a uni- 
tary transformation that diagonalizes the various correla- 
tion matrices. This transformation is called a Karhunen- 
Loeve transformation; for a stationary process with a 
correlation time much shorter than the time spanned by 
the I samples, the linear combinations of the time series 
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that diagonalize the correlation matrices asymptotically 
approach the discrete Fourier transform. Therefore we 
can approximately express our result in the frequency do- 
main where the and S^- matrices can be understood 
in terms of the power spectrum and the expectation value 
of the redshift cross-correlation respectively [cf . Eq. ( 42 ) 
and Eq. (37l]. The locally-optimal detection statistic is 
therefore 

3H'o 1 r »gw(i/i)fa'r(i/i)g:(/)g,-(/) ^^ 



Jf]„roEE''^V-2fe)5 (76) 



i — l j<i 



where f2gw(/) = £^^gw(/) E^nd fl^ — f^^a- This is the 



same optimal detection statistic S'opt of Eq. ( 66 1 (with the 
simplification of riij = 1) up to a normalization constant. 

The locally optimal statistic is optimal in the limit of 
weak signals. However, the likelihood ratio is only de- 
termined by this statistic up to a unknown factor which 
depends on the (unknown) strength of the signal. It is 
important now to distinguish between the assumed am- 
plitude of the stochastic background, e, and the true 
amplitude, etruc- The true gravitational wave spectrum 
^gw{f) is now related to the template spectrum Ug^{f) 
via rigw(/) = Cti-uc^gw- To measure the strength of the 
stochastic background given a set of pulsar observations, 
we can use the maximum likelihood estimator: the value 
of e, cmlEi for which the likelihood ratio is a maximum. 
That is, we wish to find the valu e of emle for which 
din A/c?e^|e^,LE — 0- From Eq. (74 1 we see that this esti- 
mate is 

^MLE=^-'S (77) 

where A/"^, from the O(e^) term of the log-likelihood ra- 
tio, is a normalizing factor which also includes the data. 
By substituting Eq. ( 77 1 into Eq. ( 74 1 we obtain the max- 



(78) 



imum likelihood detection statistic 
max In A 

where the terms of O(e^) have been discarded. Notice 
that this statistic is not simply the square of the cross 
correlation statistic. The data also appears in the factor 
Af~^. This factor effectively suppresses elements of the 
pulsar network where the data measured greatly exceeds 
the normal noise level. 

Some insight into the maximum likelihood detection 
statistic and the maximum likelihood amplitude estimate 
can be obtained by computing the expectation value of 
the log-likehhood ratio, Eq. (74 1. We find, to leading 
order in e. 



£-2 
true 



(S) 



V327r3 J p. 



i—l j<i 



df 



9? 

gw 



Therefore 



(l/l)fa"^r^(l/l) 

(79) 



(InA) = e^{S)-\e\N^) 



£2^.2 _ 1 2^ ( ^\ J_ 



T 



EE ^/ 



,6' 



+0{e^). (80) 



If we ignore the O(e^) terms, this is maximized when 
emle = etruc, in which case 



max (In A) 



1 



2 ^327r3 ) (3- 



I I 



EE ^/ 



1=1 j<i 



»lw(l/l)''-'"^r^(l/l) 



(81) 



This gives a scale of the value of the likelihood ratio we 
would expect to achieve. 



IV. UPPER LIMITS AND DETECTION 

Several methods exist in the LIGO literature that are 
appropriate for upper limit computation and detection 
using pulsar timing data gOl IHl HI gSJ H US] . These 
methods can be divided into two classes: Frequentist and 
Bayesian. 



We expect that the optimal statistic Eq. (66) will be 



formed from a large number of pulsar pairs. For example, 
the Parkes Pulsar Timing Array |33 ES] consists of 20 
pulsars and the optimal statistic could be constructed 
from up to 190 cross-correlation pairs. In this case we can 
make use of of the central limit theorem: The distribution 
of S'opt should be well approximated by a Gaussian with 
a mean /i = (S'opt) = ^uTq and variance cr? given by 



Eq. (68), namely. 



p(Sopt|/iCr/i) = 



exp 



-(Sopt - fJ-)"^ 



(82) 



A straightforward frequentist upper limit can then be 
set by finding the value of /iui such that in some pre- 
determined fraction C (called the confidence) of hypo- 
thetical experiments, the value of the optimal statistic 
exceeds the actual value S'opt found in the search. In 
other words we would like to find the value /Ltui such that 



dSopt p(Sopt|MulCT/l) = C. 



The solution to this is 



Sopt + \/2a^erfc-i(2(l -C)). 



(83) 



(84) 
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The assertion is that the real value of /z is less than fi^i 
with confidence C, because if /i = /Zui, a fraction C of the 
time we would have observed a value of S'opt greater than 
5'opt- An equivalent, though potentially more robust, fre- 
quentist method to set upper limits involves performing 
simulated signal injections in the timing data set. Mul- 
tiple injections are performed to determine the value of 
fiui such that a fraction C of the time the value of the op- 
timal statistic measured in the data sets with injections 
exceeds the value found in the search. Frequentist de- 
tection methods such as Neyman-Pearson or maximum- 
likchhood are well described in the literature (see, for 
example, |40j and references therein) and we will not dis- 
cuss them here. Additionally Feldman and Cousins [5T] 
provide a means to smoothly transition between upper 
limits and detection. 

Bayesian upper limits can be computed by construct- 
ing a posterior distribution using the value of the optimal 
statistic found in the search, and variance along with pri- 
ors. We begin by applying the product rule to the prob- 
ability density of fi along with the measured value S'opt 
given (T^ to write, 

pinSopt\a-f,) = p{^J.\SoptO■ii)piSopt\o■f^) 

= p{Sopt\fJ'<Jt,)p{fJ'\crii), (85) 

then solve for p(/i|S'optcr/i) to obtain Bayes' theorem. 



p(Ml'5'optO-/i) = pis opt\^J■a■^,) 



p(a^|o'a) 

P{Sopt\<Jf,) 



(86) 



the posterior probability density for /x, or equivalently 
f2a. One can then choose a prior p(/i|cr^) (for example 
requiring fj, > 0) and normalize the probability distribu- 
tion (the probability p^Soptlcp.) does not depend on fi so 
it is a prior dependent normalization constant), and find 
the /Liui such that 



(87) 



M / diip{fi\Soptcrii) = C, 



where M is the normalization constant. For sufficiently 
simple choices of the prior distribution p(^|(T^) the inte- 
gral Eq. ( 87 1 can be performed analytically to obtain the 
Bayesian analog of Eq. (84 1. As with frequentist meth- 



ods HOJ [5T] , Bayesian detection methods involve select- 
ing thresholds, in this case on the odds ratio, which is 
the ratio of the posteriors, suitably integrated over, say, 
different ranges of /i. For more details see Refs. [5^ I53j . 




Right ascension [hours] 




Right ascension [hours] 



FIG. 2: Plots of |*'-''rs=j| from Eq. (92 1 for two pulsars with 



^ = tt/2 (top panel) and ^ = tt degrees (bottom panel) 



on the sky. We begin by relaxing the assumption that 
the stochastic background is isotropic. That is, we take 



{h*AmhA'if\n')) 



327r3 



s^{n,n')SAA'S{f^f') 
xpmf\-'n,^{\f\), 



where P{Ct) is the strength or brightness [M] of gravita- 
tional waves from the direction Cl. 

In this case, the overlap reduction function takes the 
modified form. 



An 



^ / dCiP{n)F,^{n)Ff{n). (89) 



where we've ignored the pulsar phase factors, and the 
optimal filter is given by 



^'^^Qp{f) = <'^^rp»gw(|/|) 



^^?p^mpA\f\y 



(90) 



V. A PULSAR TIMING RADIOIVIETER: 
CONSTRUCTING SKYMAPS OF THE 
STOCHASTIC BACKGROUND 

A skymap may be created by computing f2gw(/) for a 
collection of pixels in the sky. We do this by assuming 
that the only signal present comes from a single location 



where we have suppressed the k index which specifies 
the particular measurement of the {ij) pulsar pair. We 
can further optimize for point sources by taking P(f2) — 
5^{Cl — Cl'). The optimal filter then becomes. 



1/1^^.(1/1)^.(1/1)^ 



(91) 
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with. 



An ^ 

A 



(92) 



Figure |2] shows two examples of the sky location depen- 
dent overlap reduction function. The top panel shows 
IFqI from Eq. (92 1 for two pulsars with ^ = tt/2 located 



at 0° Dec and 9n and 15h RA respectively. The bottom 
panel shows the same quantity for two pulsars with £, — n 
located at 0° Dec and 6h and 18h RA respectively. 

One could also imagine computing the overlap reduc- 
tion function for each term in a multipole expansion of 
P(f2). The overlap reduction function for the monopole 
term in the expansion (appropriate for an isotropic 
stochastic gravitational wave search) is t he Hellings- 
Downs curve given by Eq. ( 39 ) in Section III Surpris- 



ingly, the dipole overlap reduction function is given by a 
similarly simple equation. We find 



Tdip = (cos ai -I- cos 0:2) 
3 



cos ^ 



2 ^ / ■ C 

6 tan - In sin - 



(93) 



where as before ^ is the angle between the two pulsars, 
and ai and a2 are the angles each of two pulsars make to 
the direction of the dipole. A detailed derivation of this 
result is given in Appendix |C 2[ This result is relevant to 
searches for a dipole anisotropy in the gravitational wave 
sky using pulsar timing data. 

The sky-dependence of the sensitivity of a pulsar net- 
work can be estimated by computing the signal to noise 
for sources at the sky locations of interest. We start 
by taking the expectation value of the optimal statistic, 
Eq. ( 66 ) , using the optimal filter for a sky location f2 as- 



suming the redshift data contain a stochastic signal from 
that location. We then divide by the square root of the 
variance given in Eq. ( 68 1 . The result is proportional to 



1/2 



G(f2) 



. i—l j<i 



where we have assumed (for illustrative purposes) that 
the noise spectra of all pulsars is the same, the obser- 
vation times for all pairs is the same, and njj = 1 for 



all pulsar pairs. Figure [3] shows the Eq. (94 1, for the 
20 pulsars of the Parkes Pulsar Timing Array |3SJ 155] . 
Since most of the pulsars are in the Southern hemisphere 
the Parkes Pulsar Timing Array is most sensitive in that 
region. 

Another quantity of interest is the point spread func- 
tion, which measures the intrinsic spatial correlation of 
the skymap, or equivalently, the ability of a pulsar net- 
work to locate a stochastic source of gravitational waves. 
We construct the point spread function by computing 
the signal to noise for a source at some sky location that 
we search for using the optimal filter for some other lo- 
cation. In particular, we take the expectation value of 




Right ascension [hours] 



FIG. 3: Skymap of the sensitivity, Eq. (94 1, for the Parkes 
Pulsar Timing Array. 



the optimal statistic, Eq. ( 66 1 , using the optimal filter 



for a sky location assuming the redshift data contain 
a signal from another location fl' , then we divide by the 
square root of the variance given in Eq. ( 68 ) . The result 
is proportional to 



A{n,n') 



i—l j<.i 




(95) 



where we have again assumed that the noise spectra of all 
pulsars is the same, the observation times for all pairs is 
the same, and riij — 1 for all pulsar pairs. Figure |4| shows 



the point spread function, Eq. ( [95|, f or the 20 pulsars of 
the Parkes Pulsar Timing Arrav |351 [55] for a source at 
6h RA 45° Dec (top panel) and another at 18h RA -45° 
Dec (bottom panel). 

The point spread function can be understood in terms 
of the likelihood ratio of Sec. HID Suppose that the 



likelihood ratio is computed using the overlap reduction 
(94) function appropriate for a stochastic signal coming 
from direction when the true signal is in fact coming 
from direction Fjj. The the expectation value of the log- 
likehhood ratio is [cf. Eq. (80)] 



(InA) = e24ueEE<'^'^r,/»^-)F^,fe)C 

i—l j<i 
1=1 j<i 

+0(e') (96) 



with 



(if) (J _ 



T df 



(97) 



If ^"^i^C is approximately the same for all pulsar pairs 
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FIG. 4: Plot of the point spread function Eq. (951 for the 



Parltes pulsar timing array for a source at 6h RA 45 Dec (top 
panel) and 18h RA -45° Dec (bottom panel). 



then 



max(lnA) cx A'^{n,h'). 



(98) 



using tlie Lomb-Scargle periodogram described below, or 
modeled for each stretch. The optimally filtered data 
stretches can then be combined along the lines discussed 
in Section |IIIB[ One concern associated with breaking 
the data up into small stretches is loss of low frequency in- 
formation: Gravitational waves with periods larger than 
the length of the short stretches will be lost in this pro- 
cedure. The problem can be avoided by first computing 
the quantities Si{f)/Pi{\f\) for each of the short stretches 
and then combine them using the Dirichlet kernel to con- 
struct full time baseline versions of these quantities. 

If the spectrum is measured it can be smoothed by 
performing a running average over a small frequency win- 
dow, which if the data are stationary in the stretch the 
spectrum is estimated, is equivalent to ensemble averag- 
ing. 



B. Unevenly sampled data 

The fact that pulsar timing measurements are not 
taken continuously leads to a data set that is unevenly 
sampled in time. This poses a problem for frequency- 
domain analyses not present in their time-domain coun- 
terparts. The authors of will explore the time do- 
main approach in detail. It is unclear which of these 
two methods will turn out to be more robust. In what 
follows we address the specific issues of computing peri- 
odograms and Fourier transforms for unevenly sampled 
data sets which we think is useful in any case. 



In this sense the point spread function describes the de- 
gree to which the position of a point source of stochastic 
gravitational waves can be located in terms of the likeli- 
hood ratio. 



1. The Lomb-Scargle Periodogram 



VI. 



ISSUES WITH PULSAR TIMING DATA 



We have derived the optimal statistic for detecting a 
stochastic background. In this section we would like to 
discuss some issues associated with departures from the 
idealizations made to arrive at the optimal statistic. 



A. Colored noise and non-stationarity 

In contrast to previous methods [3 H] the techniques 
presented here do not rely on the data being white. The 
power spectra i^i(|/|) in the optimal statistic account for 
colored noise. However, the methods assume the data is 
stationary. 

If the data is non-stationary over long timcscales it 
can be divided into short stationary (or almost station- 
ary) stretches and the power spectrum can be estimated 



The problem of constructing periodograms from un- 
evenly sampled data comes up in the data analysis of vari- 
able stars. It was in precisely this context that Lomb [55] 
and Scargle [56] proposed a least-squares solution to the 
problem. The basic idea is as follows: Let x(tj) be a 
time series with zero mean sampled at i = I . . . N un- 
evenly spaced times. Now fit the time series by finding 
the coefficients Omin and that minimize the square 
of the residual 



JV 



\f) = Y,MU)-acos[2nfiU-T)]-bsm[2T:fiU-T)]y, 

(99) 



i=l 



where 



, , sin(47r/tO 
tan(47r/r) = y_Lll^ 

Ei=i cos(47r/i,) 



(100) 
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Then the periodogram is defined up to normahzation by 
the difference 



N 



1=1 

(EjIi AU) cos [27r/(t, 



(EL^(^0sin[27r/(t,-r)])^ 



(101) 



where r^in(/) is the quantity in Eq. (99 1 with a = ttmin 
and h = &min- After normahzation |57l and generalization 
to data with nonzero mean, we have the Lomb-Scargle 
periodogram 



(Eili - Mx] cos [27r/(ii - t)] j (Eili [^{ii) - Ma:] sin [27r/(ij - r) 



2a2 



Eticos^ [2^/(i.-r)] 



(102) 



where /ia; and cr^ are the mean and variance, respectively 
of x{ti). Note that the definition of r in Eq. ( 100 1 ensures 



that the resulting periodogram is independent of where 
t = 0. 



2. Fourier transforms 

The idea of using a least-squares minimization is also 
useful for constructing Fourier transforms. To do so, we 
borrow an idea from the radar community [58] . Sup- 
pose we have a timeseries, x{ti), non- uniformly sampled 
at times . . .tjq and we wish to construct its Fourier 
transform, 



i(/m) 



N 

E 

J=0 



x{tj)e 



(103) 



over M evenly spaced frequencies, /o ■ . . /m- The strat- 
egy we will employ is to use a least-squares procedure to 
find the best fit to the original timeseries after an inverse 
Fourier transform. That is, the (squared) residual to be 
minimized is given by 



N 

E 

]=0 



M 

E 

A;=0 



e(A 



(104) 



where the ^(/fe) are to be determined. Defining 

Akj=e'^''f'''\ (105) 

we can write 



= X- 



(106) 



The least-squares solution to this problem is given by 



(107) 



The problem is then purely a computational one, which, 
because of the limited amount of pulsar timing data avail- 
able, is completely tractable on a modern computer, re- 
gardless of the efficiency of the algorithm. On a final 
note, one can actually improve upon this procedure |58j 
by weighing the residual in Eq. ( |104[ ) by the square root 
of the variance, <Ta;(t^), associated with each data point 



^ 1 



^^(h) 



M 



fktj 



k=0 



(108) 



which has the advantage of automatically including the 
the error bars associated with individual pulsar timing 
data points. 



VII. SUMMARY AND CONCLUSIONS 

A stochastic background of gravitational waves could 
be detected via pulsar timing observations in the next 
5 to 10 years. This background may be astrophysical, 
such as that produced by supermassive black holes, or 
cosmological, such as that produced by a network of cos- 
mic (super)strings. In the latter case a detection would 
open a window onto a time in the early universe prior 
to recombination and could have profound consequences. 
Leveraging techniques developed for ground-based in- 
struments such as LIGO and Virgo, in this paper we have 
shown how to optimally extract the signal produced by a 
stochastic background of gravitational waves using cross- 
correlations of timing data from a pulsar timing array. 
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We started by considering the redshift induced by a 
gravitational wave on the frequency of arrival of radio 
pulses from a pulsar first derived by Detweiler [3]. The 
redshift is proportional to the difference in the metric 
perturbation at the pulsar (when a pulse is emitted) and 
at the Earth (when that pulse is received) . Using a con- 
venient coordinate independent description of the signal 
we examined the form of the signal in the frequency do- 
main. The term involving the metric perturbation at the 
pulsar is typically neglected because it can be treated as 
a sort of noise term which averages to zero in correla- 
tions of timing measurements of different pulsars. In the 
frequency domain the dependence on the metric pertur- 
bation at the pulsar is in a phase factor that depends on 
the distance to the pulsar. It is possible that if we could 
determine the distance to pulsars with sufficient accuracy 
we could use the metric perturbation at the pulsar to im- 
prove the sensitivity of our searches. Unfortunately, ac- 
curate measurements of pulsar distances are unavailable. 
By first finding the optimal cross-correlation filter, we 
have shown that for pulsar distances and gravitational 
wave frequencies typical of pulsar timing experiments, 
the metric perturbation at the pulsar can be neglected 
without a significant deviation from optimality. It is un- 
clear whether this is true for other types of gravitational 
wave searches. We have also determined the optimal way 
to combine pulsar timing data from a pulsar timing ar- 
ray, which is constructed from pairs of optimally filtered 
cross-correlations . 

We have discussed and illustrated frequentist and 
Bayesian methods for setting upper limits using the dis- 
tribution of the optimal statistic. We have shown how 
to construct a pulsar timing radiometer: A map of the 
sky created by optimizing the cross-correlation statistic 
for particular sky directions. We have also shown how to 
determine the intrinsic spatial correlation of such maps, 
which in turn determines the ability of a pulsar timing 
array to locate a source of stochastic gravitational waves. 

We have ended with a discussion of some problems 
related to realistic analysis of pulsar timing data, par- 
ticularly the issues of non-stationarity and uneven sam- 
pling. The optimal filter is constructed from power spec- 
tra of the pulsar timing data, which can be modeled or 
measured, and accounts for the effects of colored noise. 
We have described a technique, the Lomb-Scargle pe- 
riodogram, for robust spectrum estimation that can be 
used to construct the optimal filter. The optimal fil- 
ter can then be applied in the frequency domain and we 
have described a procedure for taking Fourier transforms 
of unevenly sampled data that accounts for error bars in 
the individual pulsar timing data points. The optimal fil- 
ter can also be inverse Fourier transformed and applied in 
the time domain where uneven sampling is not an issue. 
Regardless of which method turns out to be more useful 
and robust for stochastic background searches, we be- 
lieve the development of Fourier techniques for unevenly 
sampled data will be beneficial. Lommen, Romano and 
Woan 146 will examine time-domain methods in detail. 
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APPENDIX A: DERIVATION OF DETWEILER'S 
FORMULA USING THE GEODESIC EQUATION 

For completeness of presentation, we include a deriva- 
tion of Detweiler's formula. We consider, as we did in 
Section |ll] the metric perturbation due to a single gravi- 
tational wave traveling in the Cl — z direction, so Eqs. ^ 
|6]) hold. Then, if a vector is null in Minkowski space, 
the corresponding null vector, cr°, in the perturbed space- 
time gab = Vab + hab IS givcu by |[51], 



1 



(Al) 



The null vector in Minkowski space that points from 
the pulsar to the solar system is = i/(l, —a, —(3, —7), 
where, as before, a, P and 7 are the direction cosines 
with the X-, y- and z-directions, respectively. The cor- 
responding perturbed vector is readily computed from 



Eq. (All 



V 



-a(l - 



-7 



(A2) 



The geodesic equation tells us that the f-componcnt of 
CT° satisfies 

da* 



dX 



= -r 



ab^ 



(A3) 



It follows from the form of the metric perturbation in 
Eq. ^ that 



r* - - 



1 




'dgbc _^ 


dgac 


dgab 


^2 








1 . 

^dab 










/o 









1 





h+ hx 







2 







f 






\o 





0^ 





(A4) 
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where the overdot denotes a derivative with respect to t. 
The geodesic equation then reads 



da* 



1 

2' 
1 

2 

-1/ 
2 



After a httle algebra Eq. (|A2|) leads to 



- ^2 / 2 



as well as 



so that 



dv 



(A5) 



(A6) 



(A7) 



(A8) 



We proceed by writing the time derivatives in Eq. (A8) 



as derivatives with respect to the affine parameter A. In 
particular, since /i+,x = ^+.x(i — z) we have 



dh. 



dz 



dX 



dt dX 



dz dX 



(A9) 



Moreover, we also have that the frequency v = dt/dX, in 
addition to dhj^^x/dz = —dhj^/dt, and dz/dX = —1/7. 
Therefore we can write Eq. ( A9 ) as 



dh 



v{l + 7) dX 



Then Eq. (A8I, the geodesic equation, becomes 



1 dv 
V dX 



\c? - 13'^ dh+ 
2 1 + 7 ^ 



a/3 dh^ 
1 + 7^' 



which we integrate to find 
1 a2 - /32 



exp 



2 1 + 7 



a/? 
1+7" 



A/i> 



(AlO) 



(All) 



(A12) 



It is worth pointing out that the direction cosines are 
functions of the affine parameter. The dependence is in 
terms of 0{h), and we have neglected this dependence 
in going from Eq. (lAUl) to Eq. (|A12|. The final result 



is obtained by expanding this expression to first order in 



vo - v{t) 



1 a' 



where Aft,j 



2 1 + 7 



1+7 



A/ix, (A13) 



t+^x — 11'+ X ^ "■+ X is the difference between 
the metric perturbation at the pulsar and the detector. 
This expression is precisely Eq. ([t]). 



APPENDIX B: LINEARITY OF THE REDSHIFT 

In this appendix we will explicitly demonstrate that 
the total redshift is the sum of the contributions from 
gravitational waves in every direction, as written in 
Eq. (18). The derivation is merely a generalization of 



the results derived in the previous appendix. We begin 
by considering a metric perturbation, hab, in a spatial 
transverse-traceless gauge comprised of the sum of metric 
perturbations, hl'^l, from gravitational waves in N differ- 
ent directions, ^u)- Namely, 



N 



(Bl) 



where t and x form a four vector, x"" , in the background 
(Minkowski) geometry. Adjusting the notation from Ap- 
pendix A, 



dx''/dX 

iy{l,-a,-P, -7) 



(B2) 



as a null vector in Minkowski space. The null geodesic is 
perturbed by Eq. (Bl I, resulting in a 



(5s^ 



As before, our interest is in the quantity 



da' 
~dX 



-^ab<^ <^ ■ 



(B3) 



(B4) 



The spatial nature of the gauge we've chosen ensures that 



^ afc — 



1 



:9ab 



2- 



(B5) 
(B6) 



which is evident from the first line of Eq. ( A4 1 . It follows 
from these definitions that 



da* 

~dX 



= —^riabO^ o- 

= -^/iafc(.s" + <5s°)(.s^ + (5s^) 

1 ; a b 

= -^riabs s 



1 ,. 

'2 



^%(;/Vp''), 



(B7) 



where i and j are spatial indices. As before, we want 
to write the expression above in terms of the affine pa- 
rameter along s". We begin by noting that for term in 
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Eq. (|BT| 

dX 



dh^{t~n^-x) dt 



dt 



dX 



9hab{t~n(^,) ■x)d{n(^i) -x) 



dX 



dt 



dt 



V{1 + Ul^i) -p) 



dx 
dX 



dh[2{t-n^,yx) 



dt 



(B8) 



where we have used the (t — f2(j') • x) dependence of the 
metric perturbation to write the spatial derivatives as 



time derivatives along with Eq. ( B2 1 . Putting this to- 



gether with Eq. ( B7 1 , the result is that 

:h§^-n(,)-x), (B9) 



1 diy 



N 



.dX Z.2i + A(,.^- 



which can be integrated to give the redshift 



i/Q - v{t) 



N 



Y — 



(BIO) 



which is the discrete version of Eq. ( 18 ) 



APPENDIX C: THE OVERLAP REDUCTION 
FUNCTION IN THE HIGH-FREQUENCY LIMIT 



the overlap reduction function, Eq, 
the exponential factors. Thus we wis 



( 38 1 , and we ignore 
1 to evaluate 



A=+,x 



dnF^{h)F^{h) 



(Cl) 



and using the definition of F'^(f2) given by Eq. (17) we 
find 



A=+,y. ^ 



dn 



1 + • Pl 1 + • P2 



(C2) 

The two unit vectors p\ and p2 are those pointing from 
the Earth toward the first and second pulsar respec- 
tively and the polarization tensors e,J (f2) and e!^-(f2) for 

a gravitational wave traveling in direction f2 are given by 
Eqs. (2a I and (2b I respectively. To evaluate the integral 
we choose a coordinate system in which pi is parallel to 
the z-axis and P2 is in the x-z plane. Then 



Pi = (0,0,1) 

V2 = (sin ,^,0, COS 



(C3a) 
(C3b) 



1. Derivation of the Hellings-Downs curve 

In this section we derive the Hellings and Downs curve 
given by Eq. (39 1 |49]. We begin with the definition of 



where ^ is the angular separation between the two pul- 
sars. Because we have chosen coordinates in which 
p ■ rh — [cf. Eq. (3b)], the x -polarization terms vanish 



and our expression for Tq becomes 



r ^ f dCl ^ (si^^ ^ sin^ (j) — sin^ ^ cos^ cos^ — cos^ ^ sin^ 6 + 2 sin ^ cos sin 6 cos 6 cos < 
" 4 Jg2 (1 -I- cos0)(l -I- cos^cos6' -t- sin^sin^cos 0) 



(C4) 



Straightforward manipulation shows that this integral with 
becomes 



1= df2 (1 — cos 6')(1 — cosi^cos — sin^sin ^cos ( 



To =-/?(/+ J) 



(C5) 



47r ( 1 -f - cos ^ 



(C6) 
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and 



The expression in braces achieves a maximum value of 
unity when ^ = 0, so the correct normalization con- 



/■TT ^ — - ^^^^ ^.^^ 

J — — 2sin'^^ / sin6'(l — cos6')A' (C7) stant is f3 — 3/47r. With this normahzation we recover 
h Eq. (fMl). 



where we have defined 



K = 



sm 







1 + cos ^ cos 9 + sin ^ sin f? cos ( 



(C8) 



K may be trivially evaluated by contour integration in 
the complex plane. The result is 



K = 2t: 



1 + cos ^ COS 9 + I cos ^ + cos ( 



= 2tt 



sir? ^ sin^ 9 
1 =F cos C \ / 1 T cos ( 
sin^^ / V sin^e* 



2. Generalization to a Dipole Stochastic 
Background 

We now generalize the Hellings-Downs curve to the 
case of a stochastic background with a dipole moment in 
the direction D. We will start by defining the following 
quantities: 



(C9) 



where the negative sign applies when < < tt — ^ and 
the positive sign applies when tt — ^ < < tt. Hence we 
find that 



D = (sin ai cos rj, sin ai sin 77, cos ai) 

D ■ ft = cos X 

— cos ai cos 9 + sin ai sin 9 cos(0 — r/) 



(C12) 
(C13) 



J = - 47r(l - cosO 
- 47r(l + cos^) 



d9 sin ( 



(CIO) 



(C14a) 
(C14b) 



- 167r(l — cos ^) In ( sin ^ 



D ■ pi = cos a\ 

D ■ p2 = cos a2 

— cos ai cos ^ + sin ai sin ^ cos rj 

This derivation differs from the derivation of the Hellings- 
Downs curve only in that a factor D ■ Cl must be included 
in the integral. 



Combining Eqs. (|C5|, (|C6|), and ( |C10[ ), we obtain 
47r 



ro = y/5<!l + 3(l-cosO 



in I sm 2 



An 



In 
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1 — cos ^ ^ 1 
6 



DVL. 

(C15) 



rdip — 



1 



pi cos X 



This integral can be written as 

^ 

sin^ 9 (siv? ^ sin^ (j> — siv? ^ cos^ 9 cos^ 4> — cos^ ^ sin^ 9 + 2 sin ^ cos ^ sin 9 cos 9 cos < 



As in the previous section, we write 



(1 + cos 0){l + cos ^ cos 9 + sin ^ sin 6 cos ( 



Tdip- ^/3(/+ J) 



where the first term is now given by 



dCl cos x(l ~ cos ^) (1 — cos ^ cos 6 — sin ^ sin 6 cos t 



An 

= — (cos tti + cos Q!2) 

o 



(C16) 



(C17) 



(C18) 



and J is as in Eq. (C7), but K is now given by 



K = 



sin ^ cos ai cos 9 sin + (cos a2 — cos ai cos ^) sin 9 (cos (f> — sin (j)) sin 



sin^(l + cos^cos^) + sin^ ,f sin^cosf 

and may be evaluated by the same methods. The result is 



K ^ 



2tt 
sin^ ^ 



(a± cot 6* CSC 6* + b± cot^ 9 + c± csc^ 9) 



(C19) 



(C20) 
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where the following constant terms have been defined: 

cosa2 -cosaicos^ 2 ^ 

a± = cosai(l =F cos 4) ± ^ (1 ^ cos^j (L21a) 

sin ^ 

cos ao — cos «! cos f ^,9 N 

6± = Tcosai(lTcosC) ^ ^(iTcosC) (C21b) 

2 sin ^ 

cos a2 - cos ai cos ^ 2 ^ 

c± = . ^ (lip cos (C21c) 

2 sm 4 

and a+, &+, and c+ are to be used in the case where the inequality < < tt — ^ holds, and a_, and c_ are to 
be used otherwise. Thus, the integral J must again be split into two sections, and the result of the integration is 



J — An (cos ^ — 1) (cos ai + cos 0:2) — IBtt (cos ai + cos 02) tan - In |^sin - J (C22) 
Thus, we see that 

Tdip = 7r/3 (cos ai + cos 0:2) ^cos f ~ ^ ~ 4tan^ ^ In ^sin ^ (C23) 

Because we wish for Fjip to have maximal value of unity at 4 = and S, = t: (where it is clear that ai = a2 = 0), we 
must select a normalization constant of P — — 3/27r. 
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